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LONG-TERM  GOALS 

The  goal  of  this  research  is  the  development  of  a  theoretical  quantitative  understanding  of  the  dynamics  of  high  Reynolds 
number  free-surface  flows  under  large  free-surface  deformations  in  the  presence  of  one  or  multiple  surfactants.  The  long¬ 
term  goal  of  this  research  is  to  understand  the  effect  of  the  surfactants  on  (1)  turbulence,  (2)  waves,  (3)  slick  formation  and 
(4)  mass  transfer  through  the  free  surface. 

OBJECTIVES 

The  objectives  of  the  proposed  research  are  several.  First,  we  are  aiming  to  use  a  newly  developed  algorithm  for  direct 
numerical  simulation  of  free  surface  flows  with  surfactants.  This  scheme  is  capable  of  simulating  two-dimensional  flows 
with  high  Reynolds  numbers,  accommodating  nonlinear  free  surface  deformation  in  the  presence  of  one  or  multiple 
interacting  surfactants.  We  want  to  focus  on  improving  the  understanding  of  the  nonlinear  flow/free- surface/surfactant 
interactions  of  a  surfactant-contaminated  free-surface  with  a  vortex,  with  or  without  the  presence  of  waves.  Second,  we 
want  to  use  the  simulation  results  to  evaluate  a  variety  of  models  for  the  surfactant  behavior.  In  particular,  through 
comparison  with  experimental  data,  we  want  to  investigate  the  extent  of  microscopic  detail  needed  to  be  present  in  the 
surfactant  model  in  order  to  enable  an  at  least  qualitative  description  of  the  dynamic  surfactant  behavior. 

APPROACH 

The  approach  followed  for  the  accurate  and  computationally  efficient  evaluation  of  the  flow,  free-surface  deformation  and 
surfactant  concentration  involves  to  development  of  a  direct  numerical  simulation  capability  based  on  a  fully  spectral  spatial 
discretization  and  a  conformal  mapping  of  the  flow  domain  into  a  regular  parallelepiped.  The  conformal  mapping  is  key  in 
order  to  preserve  the  computational  efficiency  of  the  proposed  numerical  technique  since  it  allows  the  use  of 
computationally  efficient  fast  Poisson  algorithms  for  the  solution  of  the  Poisson  and  Helmholtz  problems  generated  from  the 
time- integration  of  the  flow  and  surfactant  concentration  evolution  equations.  The  numerical  time-integration  approach  for 
nonlinear  free-surface  flows  in  the  presence  of  insoluble  surfactants  is  a  recently  developed,  fully  implicit  scheme  based  on 
an  iterative  predictor-corrector  method.  This  approach  has  also  led  to  very  efficient  numerical  schemes  for  stationary 
boundary  problems,  where  it  has  been  extended  to  three  dimensions.  In  addition,  the  capability  of  simulating  soluble 
surfactants  can  now  be  incorporated  in  a  straightforward  manner.  We  have  also  gained  extensive  experience  in  parallel 
computing  and  have  the  capability  of  performing  our  simulations  in  powerful  parallel  computing  environments  through  very 
efficient  implementations  of  our  algorithms,  which  are  linearly  scalable  with  the  number  of  processors.  Thus,  we  are  in  a 
unique  position  to  perform  efficient  simulations  of  free-surface  flows  with  surfactants,  accommodating  fully  nonlinear  free- 
surface  deformations,  with  very  high  accuracy. 

WORK  COMPLETED 

During  the  past  year,  we  completed  the  development  of  an  algorithm  for  simulation  of  fully  non-linear  free-surface  flows  in 
the  presence  of  surfactants.  This  algorithm  exhibits  almost  linear  scalability  in  computational  cost  (Nlog2N)  with  the  number 
of  unknowns,  N.  It  is  based  on  a  fully  implicit  iterative  predictor-corrector  time-integration  scheme,  pseudoconformal 
mapping  of  the  flow-domain  to  the  computational  domain  and  a  spectrally  preconditioned  iterative  conjugate  gradient 
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solver.  Implicit  time-integration  was  necessary  for  stability  given  the  highly  nonlinear  boundary  conditions. 


The  concept  behind  the  algorithm  is  simple  and  its  implementation  is  as  follows.  An  Adams-Bashforth  scheme  is  used  for 
the  non-linear  terms  in  the  predictor  step  and  an  Adams-Moulton  scheme  for  the  corrector  step.  One  recalculates  the  non¬ 
linear  terms  at  the  current  time-step  by  repeating  the  corrector  step  in  an  iterative  fashion  till  convergence  with  high 
tolerance  (e.g.  10'12).  This  algorithm  is  only  2-3  times  more  expensive  than  the  mixed  scheme  and  has  the  advantage  that  the 
initial  guess  need  not  satisfy  the  continuity  equation  very  well,  allowing  for  more  freedom  than  a  conventional  Newton- 
Raphson  approach.  Of  course,  the  main  disadvantage  compared  to  a  full  Adams-Moulton  scheme  using  an  external 
Newton-Raphson  iteration  is  that  it  is  not  as  stable.  However,  it  is  adequately  stable  for  use  even  in  free-surface  problems. 
In  fact,  the  corrector  step  converges  faster  than  the  predictor  step,  since  the  iterative  algorithm  starts  with  a  better  initial 
guess.  Therefore,  the  algorithm  is  actually  only  a  factor  of  2  slower  than  that  using  a  mixed  time-integration  scheme. 

The  kinematic  condition  is  implemented  along  with  the  orthogonality  constraint  in  the  algorithm  that 
computes  the  time-derivatives  of  the  mapping.  This  is  achieved  by  solving  iteratively  two  separable 
Poisson  equations,  using  a  similar  method  to  that  in  Dimitropoulos  et  al.  (1998),  with  the  following 
(coupled)  boundary  conditions: 
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Once  the  time-derivatives  of  the  mapping  are  known,  we  can  then  find  the  actual  mesh  by  integrating 
in  time.  To  enforce  orthogonality  over  the  time-integration  error,  we  fit  the  top  boundary  calculated 
through  the  time-integration  to  a  Fourier  series  and  then  use  this  Fourier  series  as  a  functional 
representation  in  the  mapping  algorithm  for  a  given  surface  (described  in  Dimitropoulos  et  al,  1998). 
We  then  proceed  to  solve  for  the  flow  problem  with  a  known  mesh  at  the  new  time-step.  This  is  made 
robust  through  a  predictor-corrector  scheme.  Initially,  we  use  explicit  integration  for  the  mesh.  We 
solve  for  the  velocities  and  then  recalculate  the  time-derivatives  of  the  mesh  with  this  new  velocity  and 
use  implicit  integration  for  the  mesh.  We  proceed  to  solve  again  for  the  velocity  field.  Once  this  first 
corrector  step  is  completed,  we  repeat  it  until  the  values  we  obtain  do  not  change.  This  is  necessary 
since  we  have  decoupled  the  solution  of  the  mesh  from  the  solution  of  the  flow  problem.  This  would 
not  be  possible  using  the  fully  implicit  scheme  with  an  external  Newton-Raphson  iteration.  This 
approach  allows  us  to  solve  the  surface  diffusion  equation  for  the  surfactant  species  and  to  implement 
the  full  nonlinear  free-surface  boundary  conditions  without  needing  to  change  substantially  the 
preconditioned  conjugate  gradient  solver. 

The  surface  diffusion  equation  can  be  solved  during  each  iteration  in  a  straightforward  manner,  right 
after  the  coordinates  are  updated  through  the  mapping  algorithm.  The  equation  governing  the  transport 
of  an  insoluble  surfactant  on  an  interface  takes  the  following  form  when  expressed  in  the  orthogonal 
curvilinear  coordinate  system  that  is  utilized  in  this  work,  where  T  is  the  concentration  of  the  surfactant 
on  the  interface  and  Pes  the  surface  Peclet  number. 
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The  above  equation  is  solved  assuming  the  velocity  field  is  known,  just  as  is  performed  for  the 
mapping.  In  fact,  given  the  structure  of  the  algorithm  this  is  the  most  natural  way  to  treat  the 
convection-diffusion  problem  on  the  interface;  the  surface  concentration  is  a  property  of  the  interface, 
exactly  like  its  shape.  Finally,  due  to  the  lower  dimensionality  of  this  equation  compared  to  the 


momentum  equations  in  the  bulk,  it  is  possible  to  use  a  simple  direct  solver  and  a  fully  implicit  Adams- 
Moulton  time-integration  scheme  to  solve  the  system  of  algebraic  equation  that  results  after 
pseudospectral  discretization  using  Fourier  collocation.  The  cost  of  solving  this  equation  does  not 
affect  the  overall  computational  cost  of  the  algorithm. 

Once  the  surface  concentration  is  known,  one  can  calculate  the  surface  tension  through  a  constitutive 
relationship.  We  have  chosen  to  use  a  nonlinear,  Langmuir  type  relationship  which  is  adequate  to 
capture  local  surfactant  saturation  effects.  The  equation  considered  is 

(7  =  1  +  Ma{\. + y3)ln 


where  Ma  is  the  Marangoni  number,  which  represents  the  ratio  of  the  interfacial  tension  gradient  forces 
to  bulk  viscosity  forces  and  the  (3  is  a  parameter  that  represents  nonlinear  packing  effects.  When  p  is 
very  large,  ideal  behavior  (linear  variation  of  surface  forces  with  concentration)  is  recovered. 


For  the  free-surface  boundary  conditions,  we  utilize  the  Boussinesq-Scriven  model  for  the  surface 
stress,  which  assumes  linear  dependence  of  the  stress  on  the  strain  rate.  When  we  express  the  boundary 
conditions  in  the  pseudoconformal  coordinate  system  utilized  in  this  work  they  take  the  form: 
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tangential  stress, 


where  Ca  is  the  capillary  number,  Bo  is  the  Boussinesq  number  and  H  is  the  mean  curvature  of  the 
interface.  We  implement  these  boundary  conditions  in  the  preconditioner,  along  with  the  divergence- 
free  condition  after  transforming  them  into  the  following  set  of  equations: 
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The  constant  K  is  the  average  of  +  HBo)>  RBC1,  RBC2,  RBC3  are  the  residuals  of  the  boundary  conditions  calculated 

in  the  conjugate  gradient  level  and  i  is  the  index  of  the  Concus  and  Golub  type  iteration  performed  in  the  preconditioner. 
This  scheme  is  constructed  so  that  the  influence  matrix  method  can  be  efficiently  implemented  in  order  to  enforce  the 
divergence-free  condition. 


RESULTS 


We  have  validated  the  algorithm  by  comparing  its  predictions  against  those  of  linear  theory  of  wave  damping  in  the 
presence  of  surfactants.  This  work  will  be  communicated  in  a  series  of  forthcoming  articles  (see  Dimitropoulos  1999c) 
which  will  include  the  methodology  and  results  in  a  comprehensive  fashion.  Here,  we  present  briefly  results  from 
simulations  of  a  purely  elastic  interface  without  viscous  or  molecular  diffusion  effects.  In  the  case  that  surfactants  are 
present,  two  modes  of  motion  exist,  transverse  as  well  as  longitudinal  waves.  Tables  1  and  2  present  a  comparison  of  linear 
theory  with  simulation  results  for  damping  of  transverse  and  longitudinal  waves.  The  initial  amplitude  was  equal  to  10'3  in  a 
channel  2  units  deep  and  1  unit  wide,  Re=104  (based  on  the  half-depth  and  the  gravity),  Ca=10'1  (corresponding  to  water), 
Ma=4/9,  Bo=0,  Pe  =1018  and  the  time-step  was  equal  to  10'3. 


Linear  Theory 

Numerical  Simulation 

Damping  Factor 

-0.062571 

-0.062599 

Period  of  Oscillation 

2.43798 

2.43815 

Table  1:  Damping  of  transverse  waves:  Comparison  of  simulations  with  linear  theory. 


Linear  Theory 

Numerical  Simulation 

Damping  Factor 

-1.222236 

-1.298954 

Period  of  Oscillation 

2.98021 

2.88204 

Table  2:  Damping  of  longitudinal  waves:  Comparison  of  simulations  with  linear  theory. 

It  can  be  readily  observed  that  the  simulations  agree  well  with  theory.  It  should  be  noted  that  the  larger  discrepancy  in  the 
results  for  longitudinal  waves  is  due  to  the  sensitivity  of  the  simulations  in  the  initial  guess  that  makes  the  initial  motion, 
which  is  mostly  of  longitudinal  character,  to  evolve  and  include  a  small  superimposed  contribution  of  the  transverse  mode. 
We  must  emphasize  here  that  we  found  that  the  flow  field  that  corresponds  to  each  mode  is  quite  different.  The  horizontal 
velocity  for  longitudinal  waves  in  the  vicinity  of  the  surface  is  one  order  of  magnitude  greater  than  for  transverse  waves.  As 
a  result,  the  pressure  distribution  also  has  distinct  differences,  which  persist  further  away  from  the  surface  than  the  velocity 
differences.  Detailed  results  that  clearly  demonstrate  the  different  physics  for  each  type  of  motion  are  included  in  a  recent 
doctoral  dissertation  supported  by  this  project  and  will  be  presented  in  a  series  of  future  publications. 

IMPACT/APPLICATIONS 

We  have  developed  a  computationally  efficient  numerical  scheme  for  solving  spectrally,  two-dimensional,  time-dependent 
free-surface  flows  in  the  presence  of  insoluble  surfactants.  The  implementation  was  such  that  the  almost  linear  scalability 
and  exponential  convergence  with  mesh  refinement  characteristics  are  retained.  An  orthogonal  mapping  algorithm,  an 
efficient  and  robust  iterative  solver  incorporating  the  influence  matrix  method  for  satisfying  the  incompressibility  condition 
and  an  iterative  predictor-corrector  implementation  have  been  developed.  The  algorithm  is  parallelizable  in  a 
straightforward  fashion  and  exhibits  linear  scalability  of  its  performance  with  the  number  of  processors  used.  This  is  the 
first  time  in  our  knowledge  that  a  fully  spectral  method  for  simulation  of  fully  nonlinear  free-surface  flows  with  surfactants 
has  been  developed.  This  method  has  unique  capabilities  to  elucidate  the  complex  phenomena  that  occur  in  situations  where 
there  is  nonlinear  coupling  of  fluid  flow,  interface  deformation  and  mass  transport  of  surfactants. 

TRANSITIONS 

The  main  development  of  this  work  along  with  that  performed  over  the  past  5  years  is  a  quite  general  in  applicability  set  of 
algorithms,  which  extend  the  use  of  fully  spectral  methods  to  flows  within  moderately  complex  geometries.  The  tools 
developed  in  this  work  i.e.  the  pseudoconformal  mapping  and  the  efficient  spectral  preconditioner  are  also  applicable  to 
other  computational  fluid  mechanics  and  transport  phenomena  problems.  All  of  these  are  expected  to  find  significant  use 
within  the  computational  fluid  dynamics  community,  since  their  possible  applications  span  a  quite  wide  area. 

RELATED  PROJECTS 

The  core  of  the  spectral  numerical  method,  after  suitable  modification,  has  been  successfully  used  in  direct  numerical 
simulations  (DNS)  of  turbulent  viscoelastic  flows  in  a  joint  University  of  Delaware  -  NRL  effort  for  the  investigation  of 


polymer-induced  turbulent  drag  reduction.  There  has  been  substantial  progress  after  our  latest  simulations  in  the  CRAY- 
0rigin2000  at  National  Center  for  Supercomputing  Applications  (NCSA).  This  work  has  produced  so  far  four  publications 
(Sureshkumar  et  al .,  1997;  Dimitropoulos  et  al. ,  1998b,  Beds  and  Dimitropoulos,  1999a,  Dimitropoulos  et  al ,  1999b)  and  is 
being  continued  to  study  further  the  phenomenon  of  drag  reduction.  The  implementation  in  multiprocessor  environments  of 
the  codes  developed  in  the  current  project  has  been  greatly  facilitated  from  the  above-mentioned  experience. 
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